

#include <algorithm>
#include <chrono>
#include "femSolids.h"
#include "log.h"
#include "util.h"
#include <boost/algorithm/string.hpp>
#include <string>
#include <unordered_map>
#include <unistd.h>
#include "MyTime.h"
#include "TimeFunction.h"
#include "utilitiesmaterial.h"
#include "NewMaterial.h"
#include "utilitiesmaterial.h"
#include "NewLagrangeElement.h"
#include "QuadratureUtil.h"
#include "BasisFunctionsLagrange.h"



using namespace std;

extern  std::vector<unique_ptr<TimeFunction> > timeFunction;
extern  MyTime           myTime;
extern  bool  debug;

femSolids* globalFemSolidsInstance = nullptr;
femSolids::femSolids(string& file_name)
{
    input_file_name = file_name;
    globalFemSolidsInstance = this;
    MPI_Comm_size(MPI_COMM_WORLD, &n_mpi_procs);
    MPI_Comm_rank(MPI_COMM_WORLD, &this_mpi_proc);

    nNode_owned  = 0;
    nElem_global = 0;
    nNode_global = 0;

    nNode_local = 0;    
    fileCount = 0;
    ndof = 0;

    numElemTypes = 0;
    numMaterials = 0;

    dispDOF = 0;
    presDOF = 0;
    totalDOF = 0;

    intVarFlag = false;
    IMPLICIT_SOLVER = true;
    MIXED_ELEMENT = false;
    convergedFlagPrev = convergedFlag = false;


    loadStepConverged = 0;
    loadFactor = 0.0;
    loadFactorPrev = 0.0;
    loadFactorPrev2 = 0.0;
    loadfactorVec.push_back(loadFactor);

    localStiffnessError = 0;

    NONLINEAR_SOLVER_TYPE = SOLVER_TYPE_NEWTONRAPHSON;

    AlgoType = 2;

    conv_tol = 1.0e-6;

    elems = nullptr;

    dirname = get_current_dir_name();
    vector<string> stringlist;
    boost::algorithm::split(stringlist, dirname, boost::is_any_of("/"), boost::token_compress_on);
    for(auto& str: stringlist)  boost::trim(str);

    dirname = stringlist[stringlist.size()-1];
}


femSolids::~femSolids()
{
  if(elems != NULL)
  {
    for(int ii=0;ii<nElem_local;ii++)
      delete elems[ii];

    delete [] elems;
    elems = NULL;
  }
}




int  femSolids::deallocatePetscObjects()
{
  solverPetsc->free();

  return 0;
}

[[nodiscard]] COMMON_ERR_CODE femSolids::read_file_op(ifstream& infile, string& op_type) {
    COMMON_ERR_CODE err_code; 
    if(op_type == "read_common") {
        err_code = readCommonInfo(infile);
    }
    else if(op_type == "read_element_pre") {
        err_code = readSolidElements(infile);
    }
    else if(op_type == "read_element_local") {
        err_code = readSolidElementsLocal(infile);
    }
    else if(op_type == "read_node_local") {
        err_code = readSolidNodesLocal(infile);
    }
    else {
        err_code = NO_READ_FILE_OP;
    }
    return err_code;
}


[[nodiscard]] COMMON_ERR_CODE femSolids::read_file_op1(ifstream& infile, string& op_type) {
    COMMON_ERR_CODE err_code; 
    if(op_type == "read_common2") {
        err_code = readCommonInfo2(infile);
    }else{
        err_code = NO_READ_FILE_OP;
    }
    return err_code;
} 


[[nodiscard]] COMMON_ERR_CODE femSolids::readCommonInfo2(ifstream& infile) {
    COMMON_ERR_CODE err_code;
    string line;
    vector<string>  stringlist;

    while(getline(infile,line))
    {
        boost::trim(line);

        if( isAValidLine(line) )
        {
            if( (line.compare(string("Nodes")) == 0) || (line.compare(string("nodes")) == 0) )
            {
                cout << "Nodes" << endl;
                readNodesLocal(infile, line);
            }
        }
    }
    err_code = COMMON_SUCCESS;

    return err_code;
}




[[nodiscard]] COMMON_ERR_CODE femSolids::readCommonInfo(ifstream& infile) {
    COMMON_ERR_CODE err_code;
    cout << " femSolids::readInput " << endl;

    string line;
    vector<string>  stringlist;

    while(getline(infile,line))
    {
        //if (line == "") continue;
        boost::trim(line);

        //std::cout << line << std::endl;

        if( isAValidLine(line) )
        {
            //std::cout << line << std::endl;

            if(line.compare(string("Model Type")) == 0)
            {
                cout << "Model Type" << endl;
                readModelType(infile, line);
            }
            else if( (line.compare(string("Nodes")) == 0) || (line.compare(string("nodes")) == 0) )
            {
                cout << "Nodes" << endl;
                readNodes(infile, line);
            }
            else if( (line.compare(string("Prescribed Boundary Conditions")) == 0) || (line.compare(string("prescribed boundary conditions")) == 0) )
            {
                cout << "Prescribed Boundary Conditions" << endl;
                readPrescribedBCs(infile, line);
            }
            else if(line.compare(string("Face Load Elements")) == 0)
            {
                cout << "Face Load Elements" << endl;
                readSurfaceElements(infile, line);
            }
            else if( (line.compare(string("Nodal Forces")) == 0) || (line.compare(string("nodal forces")) == 0) )
            {
                cout << "Nodal Forces" << endl;
                readNodalForces(infile, line);
            }
            else if(line.compare(string("Element")) == 0)
            {
                cout << "Element Type" << endl;
                readElementProps(infile, line);
            }
            else if(line.compare(string("Material")) == 0)
            {
                cout << "Material Type" << endl;
                readMaterialProps(infile, line);
            }
            else if(line.compare(string("Body Force")) == 0)
            {
                cout << "Body Force" << endl;
                readBodyForce(infile, line);
            }
            else if(line.compare(string("Time Function")) == 0)
            {
                cout << "Time Function" << endl;
                readTimeFunctions(infile, line);
            }
            else if(line.compare(string("Solver")) == 0)
            {
                cout << "Solver" << endl;
                readSolverDetails(infile, line);
            }
            else if(line.compare(string("Nodal Data Output")) == 0)
            {
                cout << "Output" << endl;
                readNodalDataForOutput(infile, line);
            }
            else if(line.compare(string("Output")) == 0)
            {
                cout << "Output" << endl;
                readOutputDetails(infile, line);

            }
        }
    }
    err_code = COMMON_SUCCESS;

    return err_code;
}


[[nodiscard]] COMMON_ERR_CODE femSolids::prepareInputData(){
    COMMON_ERR_CODE err_code;
    
    nElem_local_max = nElem_local;
    vector<int> nElem_local_count(n_mpi_procs, 0);

    /*** Send nElem_local count ***/
    for (int pid_recv = 0; pid_recv < n_mpi_procs; pid_recv++) {
        if(pid_recv != this_mpi_proc) {
            MPI_Send(&nElem_local, 1, MPI_INT, pid_recv, MPI_ELE_MAXCOUNT, MPI_COMM_WORLD);

            LOG_INFO("Send nElem_local from pid[" + to_string(this_mpi_proc) + "] to pid[" + to_string(pid_recv) +
                "]\tcount: " + to_string(nElem_local)
                + string(__FILE__) + ":" + to_string(__LINE__));
        }
    }
    /*** Recv nElem_local count ***/
    for (int pid_send = 0; pid_send < n_mpi_procs; pid_send++) {
        if(pid_send != this_mpi_proc) {
            MPI_Recv(&nElem_local_count[pid_send], 1, MPI_INT, pid_send, MPI_ELE_MAXCOUNT, MPI_COMM_WORLD, MPI_STATUS_IGNORE);

            LOG_INFO("Recv nElem_local from pid[" + to_string(pid_send) + "] to pid[" + to_string(this_mpi_proc) +
                "]\tcount: " + to_string(nElem_local_count[pid_send])
                + string(__FILE__) + ":" + to_string(__LINE__));
        }
    }
    /*** Comp nElem_local_max count ***/
    for (int proc = 0; proc < n_mpi_procs; proc++) {
        nElem_local_max = max(nElem_local_max, nElem_local_count[proc]);
    }
    LOG_INFO("Recv max end, nElem_local_max: " + to_string(nElem_local_max)
        + string(__FILE__) + ":" + to_string(__LINE__));

    FILE *file;
    char path[50];
    long pages, dummy, data, stack, text, total_vm, rss;
    sprintf(path, "/proc/%d/statm", getpid());

    int ii, ee, ind, matId, elemId, npElem;

    ind = nNode_global*ndof;

    //wj
    int ind1 = num_nodeInLocalElem * ndof;



    LOG_INFO("Start init solndata, nNode_global: " + to_string(nNode_global)
        + " ndof: " + to_string(ndof) + string(__FILE__) + ":" + to_string(__LINE__));
    SolnData.initialise(ind, nNode_global, 0, 0);

    //wj
    SolnData.initialise1(ind1, nNode_global, 0, 0);   
    LOG_INFO("init solndata end");

    file = fopen(path, "r");
    if (file != NULL) {
        fscanf(file, "%ld %ld %ld %ld %ld %ld",
           &total_vm, &data, &stack, &text, &dummy, &rss);
        fclose(file);

        pages = sysconf(_SC_PAGESIZE) / 1024;
        double rss_mb = rss * pages / 1024.0;

        LOG_INFO("prepare node data end RSS in usage: " + to_string(rss_mb * 1e-3) + " Gbytes");
    }

    vector<int>  nodeNums;

    //vector<int>  nodeInLocalElem;
    num_nodeInLocalElem2 = 0; 

    LOG_INFO("Start prepare ele data, ele size: " + to_string(nElem_local)
        + " array size: " + to_string(elemConn.size())
        + string(__FILE__) + ":" + to_string(__LINE__));

    elems = new ElementBase* [nElem_local];
    for(ee = 0; ee < nElem_local; ++ee) {
      npElem = elemConn[ee].size()-3;

      nodeNums.resize(npElem);
      nodeInLocalElem.resize(npElem);
      nNode_local = nNode_owned;

      for(ii=0; ii<npElem; ii++){

        nodeNums[ii] = elemConn[ee][3+ii];
      
        if (node_map_g2l.find(nodeNums[ii]) == node_map_g2l.end())
        {
            if(nodeNums[ii]>= node_start && nodeNums[ii] <= node_end)
            {
                nodeInLocalElem[ii] = nodeNums[ii] - node_start;
                node_map_g2l.insert({nodeNums[ii],nodeInLocalElem[ii]});
                node_map_l2g[nodeInLocalElem[ii]] = nodeNums[ii];
            }
            else
            {
                nodeInLocalElem[ii] = nNode_local + num_nodeInLocalElem2;
                node_map_g2l.insert({nodeNums[ii],nodeInLocalElem[ii]});
                node_map_l2g[nodeInLocalElem[ii]] = nodeNums[ii];
                num_nodeInLocalElem2++;
            }
        }
        else
        {
            nodeInLocalElem[ii] = node_map_g2l[nodeNums[ii]];
        }
    }
      elemId = elemConn[ee][0];
      matId  = elemConn[ee][1];

      elems[ee] = NewLagrangeElement(ElementTypeDataList[elemId]->getElemTypeNameNum());

      elems[ee]->elenum       = elemConn_idx[ee];
      elems[ee]->nodeNums     = nodeNums;

      elems[ee]->nodeInLocalElem = nodeInLocalElem;
      elems[ee]->nodeNums     = nodeNums;     
      elems[ee]->SolnData     = &(SolnData);
      elems[ee]->GeomData     = &(GeomData);
      elems[ee]->MatlData     = MatlDataList[matId];
      elems[ee]->ElemTypeData = ElementTypeDataList[elemId];

      elems[ee]->prepareElemData();
}
// for (auto it = node_map_g2l.begin(); it != node_map_g2l.end(); ++it) {

//     LOG_INFO("Global Node: " + to_string(it->first));
//     LOG_INFO("Local Node: " + to_string(it->second));

// }


    LOG_INFO("Prepare ele data end " + string(__FILE__) + ":" + to_string(__LINE__));

    PetscPrintf(MPI_COMM_WORLD, "\n elements are created and prepared \n");

    err_code = COMMON_SUCCESS;

    PetscPrintf(MPI_COMM_WORLD, "     femSolids::prepareInputData()  .... FINISHED ...\n\n");

    return err_code;
}




int femSolids::readConfiguration(string& fname)
{
    cout << " femSolids::readConfiguration " << endl;

    ifstream  infile(fname);

    if(infile.fail())
    {
       cout << " Could not open input file " << endl;
       exit(-1);
    }


    string line;
    vector<string>  stringlist;

    while(getline(infile,line))
    {
        boost::trim(line);
 
        //cout << " size       = " << line.size() << endl;
        //cout << " first word = " << line[0] << endl;

        if( isAValidLine(line) )
        {
            //std::cout << line << std::endl;

            if(line.compare(string("files")) == 0)
            {
                cout << "files" << endl;
                // read the files
                //readInputFiles(infile, line);
            }
            else if(line.compare(string("Material")) == 0)
            {
                cout << "Material" << endl;
                readMaterialProps(infile, line);
            }
            else if(line.compare(string("Body Force")) == 0)
            {
                cout << "Body Force" << endl;
                readBodyForce(infile, line);
            }
            else if(line.compare(string("Time Functions")) == 0)
            {
                cout << "Time Functions" << endl;
                readTimeFunctions(infile, line);
            }
            else if(line.compare(string("Solver")) == 0)
            {
                cout << "Solver" << endl;
                readSolverDetails(infile, line);
            }
            else if(line.compare(string("Initial Conditions")) == 0)
            {
                cout << "Initial Conditions" << endl;
                readInitialConditions(infile, line);
            }
            else if(line.compare(string("Nodal Forces")) == 0)
            {
                cout << "Nodal Forces" << endl;
                readNodalForces(infile, line);
            }
            else if(line.compare(string("Output")) == 0)
            {
                cout << "Output" << endl;
                //readOutputDetails(infile, line);
            }
            else
            {
                cout << "key =  " <<  line << endl;
                throw runtime_error("Key not found in femSolids::readConfiguration ...");
                //return -1;
            }
      }
    }

    infile.close();

    cout << " Configuration file is successfully read " << endl;

    return 0;
}





int  femSolids::readModelType(ifstream& infile, string& line)
{
    vector<string>  stringlist;

    // read {
    getline(infile,line);    boost::trim(line);

    cout << line << endl;

    while( infile && (line != "}") )
    {
        getline(infile,line);    boost::trim(line);

        if(line[0] == '}') break;


        if( isAValidLine(line) )
        {
            boost::algorithm::split(stringlist, line, boost::is_any_of(":"), boost::token_compress_on);
            for(auto& str: stringlist)  boost::trim(str);

            if(stringlist[0] == "dimension")
            {
                ndim = stoi(stringlist[1]);

                assert( (ndim > 0) && (ndim < 4) );

                ndof = ndim;
            }
            else if(stringlist[0] == "behaviour")
            {
                if( (stringlist[1] == "pstress") || (stringlist[1] == "PSTRESS") )
                  GeomData.setModelBehaviour(PSTRESS);
                else if ( (stringlist[1] == "pstrain") || (stringlist[1] == "PSTRAIN") )
                  GeomData.setModelBehaviour(PSTRAIN);
                else if ( (stringlist[1] == "axisym") || (stringlist[1] == "AXISYM") )
                  GeomData.setModelBehaviour(AXISYM);
                else
                {
                  throw runtime_error("Behaviour type not available in femSolids::readModelType");
                }
            }
            else if(stringlist[0] == "thickness")
            {
                GeomData.setThicnkess(stod(stringlist[1]));
            }
            else
            {
                throw runtime_error("Option not available in femSolids::readModelType");
            }
        }
    }

    return 0;
}



int  femSolids::readNodesLocal(ifstream& infile, string& line)
{
    vector<string>  stringlist;

    // read the number of nodes
    getline(infile,line);    boost::trim(line);

    nNode_global = stoi(line);


    vector<myPoint>  nodePosDataLocal;
    nodePosDataLocal.resize(num_nodeInLocalElem);

    int m = 0;


    
    for(int ii=0; ii<nNode_global; ii++)
    {
        getline(infile,line);        boost::trim(line);

        boost::algorithm::split(stringlist, line, boost::is_any_of("\t "), boost::token_compress_on);
        for(auto& str: stringlist)  boost::trim(str);

        // LOG_INFO("positions: "      + to_string(stod(stringlist[1])) + ":"+ to_string(stod(stringlist[2])) +":"+
        // string(__FILE__) + ":" + to_string(__LINE__));
        

        if(node_map_g2l.find(node_map_get_new[stoi(stringlist[0])-1]) != node_map_g2l.end())
        {
            // LOG_INFO("stringlist: "      + to_string(node_map_get_new[stod(stringlist[0])-1]) +  
            // string(__FILE__) + ":" + to_string(__LINE__));

            m++;
            int local_pos_idx = node_map_g2l[node_map_get_new[stod(stringlist[0])-1]];
            nodePosDataLocal[local_pos_idx] [0] = stod(stringlist[1]); 
            nodePosDataLocal[local_pos_idx] [1] = stod(stringlist[2]);           
            if(ndim == 3) nodePosDataLocal[local_pos_idx][2] = stod(stringlist[3]);
            // LOG_INFO("nodePosDataLocal: "      + to_string(nodePosDataLocal[local_pos_idx] [0]) +  to_string(nodePosDataLocal[local_pos_idx] [1])  + 
            // string(__FILE__) + ":" + to_string(__LINE__));
        }
    }



    GeomData.setDimension(ndim);

    GeomData.setNodalPositionsLocal(nodePosDataLocal);

  
    myPoint  pt;

    for(int ii=0;ii<num_nodeInLocalElem;ii++)
    {

        pt =  GeomData.NodePosOrigLocal[ii];


        //std::cout << "GeomData.NodePosOrigLocal[" << ii << "] = " << GeomData.NodePosOrigLocal[ii] << std::endl; 
        GeomData.NodePosCurLocal[ii]  = pt;
        GeomData.NodePosNewLocal[ii]  = pt;

    }

    //GeomData.updateNodesAfterPartition();
    return 0;
}



int  femSolids::readNodes(ifstream& infile, string& line)
{
    vector<string>  stringlist;

    // read the number of nodes
    getline(infile,line);    boost::trim(line);

    nNode_global = stoi(line);

    cout << "nNode_global = " << nNode_global << endl;
    cout << "ndim = " << ndim << endl;

    // vector<myPoint>  nodePosData;
    // nodePosData.resize(nNode_global);

    // for(int ii=0; ii<nNode_global; ii++)
    // {
    //     getline(infile,line);        boost::trim(line);

    //     boost::algorithm::split(stringlist, line, boost::is_any_of("\t "), boost::token_compress_on);
    //     for(auto& str: stringlist)  boost::trim(str);

    //     nodePosData[ii][0] = stod(stringlist[1]);
    //     nodePosData[ii][1] = stod(stringlist[2]);
    //     if(ndim == 3) nodePosData[ii][2] = stod(stringlist[3]);
    //     // LOG_INFO("nodePosData: "      + to_string(nodePosData[ii] [0]) +  to_string(nodePosData[ii] [1])  +  to_string(ii)  + 
    //     // string(__FILE__) + ":" + to_string(__LINE__));
    // }
    // // LOG_INFO("nodePos: "      + to_string(nodePosData.size()) + 
    // // string(__FILE__) + ":" + to_string(__LINE__));
    // GeomData.setDimension(ndim);
    // GeomData.setNodalPositions(nodePosData);

    return 0;
}





[[nodiscard]] COMMON_ERR_CODE femSolids::readSolidElements(ifstream& infile) {
    COMMON_ERR_CODE err_code;
    string line;

    while(getline(infile,line)) {
        boost::trim(line);
        if(isAValidLine(line)) {
            if((line.compare(string("Elements")) == 0) || (line.compare(string("elements")) == 0)) {
                getline(infile,line);    
                boost::trim(line);

                nElem_global    = stoi(line);
                ele_local_size  = nElem_global / n_mpi_procs;
                nElem_start     = this_mpi_proc * ele_local_size;
                nElem_end       = (this_mpi_proc != (n_mpi_procs - 1)) ? (nElem_start + ele_local_size) : nElem_global; 
                nElem_local     = nElem_end - nElem_start;

                LOG_INFO("nElem_global = " + to_string(nElem_global) + string(__FILE__) + ":" + to_string(__LINE__));
                LOG_INFO("nElem_local  = " + to_string(nElem_local ) + string(__FILE__) + ":" + to_string(__LINE__));
                LOG_INFO("nElem_start  = " + to_string(nElem_start ) + string(__FILE__) + ":" + to_string(__LINE__));
                LOG_INFO("nElem_end    = " + to_string(nElem_end   ) + string(__FILE__) + ":" + to_string(__LINE__));

                if(nElem_local > ELE_LOC_LIMIT) {
                    err_code = OUT_OF_MEMORY;
                    break;
                }

                int cout_start = nElem_start;
                while(cout_start > 0) {
                    getline(infile,line);
                    cout_start--;
                }

                vector<string>  stringlist;
                elemConn.resize(nElem_local);

                for(int ee = 0; ee < nElem_local; ee++) {
                    getline(infile,line);
                    boost::trim(line);
                    boost::algorithm::split(stringlist, line, boost::is_any_of("\t "), boost::token_compress_on);

                    int ind = stringlist.size();
                    elemConn[ee].resize(ind - 1);

                    for(int j = 1; j < ind; j++) {
                        elemConn[ee][j - 1] = stoi(stringlist[j]) - 1;
                    }
                }

                err_code = COMMON_SUCCESS;
                break;
            }
        }
    }

    return err_code;
}

[[nodiscard]] COMMON_ERR_CODE femSolids::readSolidElementsLocal(ifstream& infile) {
    COMMON_ERR_CODE err_code;
    string line;

    while(getline(infile,line)) {
        boost::trim(line);
        if(isAValidLine(line)) {
            if((line.compare(string("Elements")) == 0) || (line.compare(string("elements")) == 0)) {
                getline(infile,line);    

                nElem_local = ele_local_idx_set.size();
                if(nElem_local > ELE_LOC_LIMIT) {
                    err_code = OUT_OF_MEMORY;
                    break;
                }
                LOG_INFO("Read element local start size: " + to_string(nElem_local)
                    + string(__FILE__) + ":" + to_string(__LINE__));

                int local_idx = 0;
                vector<string>  stringlist;
                elemConn.resize(nElem_local);
                elemConn_idx.resize(nElem_local);

                for(int global_idx = 0; global_idx < nElem_global; global_idx++) {
                    if(local_idx == nElem_local) {
                        err_code = COMMON_SUCCESS;
                        break;
                    }

                    getline(infile,line);
                    if(ele_local_idx_set.find(global_idx) == ele_local_idx_set.end()) {
                        continue;
                    }
                    LOG_DEBUG("pid[" + to_string(this_mpi_proc) + "]\t<-- ele[" + to_string(global_idx) + "]");

                    boost::trim(line);
                    boost::algorithm::split(stringlist, line, boost::is_any_of("\t "), boost::token_compress_on);

                    int ind = stringlist.size();
                    elemConn[local_idx].resize(ind - 1);

                    elemConn_idx[local_idx] = global_idx;
                    for(int j = 1; j < ind; j++) {
                        elemConn[local_idx][j - 1] = stoi(stringlist[j]) - 1;
                    }

                    local_idx++;
                }
                
                if(local_idx != nElem_local) {
                    err_code = READ_ELE_LESS;
                    break;
                }

                err_code = COMMON_SUCCESS;
                break;
            }
        }
    }

    return err_code;
}

// TODO::Parallel & remap of nodes
// 

[[nodiscard]] COMMON_ERR_CODE femSolids::readSolidNodesLocal(ifstream& infile) {
    COMMON_ERR_CODE err_code;
    node_map_get_old.resize(nNode_global, 0);
    node_map_get_new.resize(nNode_global, 0);
    vector<int> node_proc_id(nNode_global, -1);

    //wj
    num_nodeInLocalElem = 0;

    for(int ee = 0; ee < nElem_local; ee++) {
        int npElem = elemConn[ee].size() - 3;

        for (int i = 0; i < npElem; i++) {
            node_map_get_old[elemConn[ee][3 + i]]++;


            if(node_map_g2l.find(elemConn[ee][3 + i]) == node_map_g2l.end())
            {
                node_map_g2l.insert({elemConn[ee][3 + i], 0});
                num_nodeInLocalElem++;
            } 
        }
    }
    node_map_g2l.clear();
    node_map_l2g.resize(num_nodeInLocalElem);
    // LOG_INFO("node_map_l2g: " + to_string(num_nodeInLocalElem));
    if(this_mpi_proc == 0) {
        for(int i = 0; i < nNode_global; i++) {
            if(node_map_get_old[i] != 0) {
                node_proc_id[i] = 0;
            }
        }
    }

    MPI_Barrier(MPI_COMM_WORLD);

    if(this_mpi_proc != 0) {
        MPI_Send(&node_map_get_old[0], nNode_global, MPI_INT, 0, MPI_NODE_MAP_BG, MPI_COMM_WORLD);
    }

    if(this_mpi_proc == 0) {
        vector<int> node_map_buffer(nNode_global, 0);
        for(int pid_send = 1; pid_send < n_mpi_procs; pid_send++) {
            LOG_INFO("recv node map start from pid: " + to_string(pid_send));
            MPI_Recv(&node_map_buffer[0], nNode_global, MPI_INT, pid_send, MPI_NODE_MAP_BG, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
            for(int i = 0; i < nNode_global; i++) {
                if(node_map_buffer[i] > node_map_get_old[i]) {
                    node_map_get_old[i] = node_map_buffer[i];
                    node_proc_id[i] = pid_send;
                }
            }
            LOG_INFO("recv node map end from pid: " + to_string(pid_send));
        }
    }

    if(this_mpi_proc == 0) {
        int left, right, mid_pid;
        for(left = 0; left < nNode_global; left++) {
            if(node_proc_id[left] == -1) {
                right = left;
                while((right < nNode_global) && (node_proc_id[right] == -1)) {
                    right++;
                }
                if(left == 0) {
                    mid_pid = node_proc_id[right];
                } else if(right == nNode_global) {
                    mid_pid = node_proc_id[left - 1];
                } else {
                    mid_pid = (node_proc_id[left - 1] + node_proc_id[right]) / 2;
                }
                for(int i = left; i < right; i++) {
                    node_proc_id[i] = mid_pid;
                }
                left = right;
            }
        }

        for(int i = 0; i < nNode_global; i++) {
            LOG_DEBUG("node[" + to_string(i) + "]\t--> " + to_string(node_proc_id[i]));
        }
    }

    MPI_Barrier(MPI_COMM_WORLD);

    MPI_Bcast(&node_proc_id[0], nNode_global, MPI_INT, 0, MPI_COMM_WORLD);

    nNode_owned = std::count(node_proc_id.begin(), node_proc_id.end(), this_mpi_proc);

    std::vector<int>  nodelist_owned(nNode_owned);

    int kk = 0;
    for(int ii = 0; ii < nNode_global; ii++)
    {
        if(node_proc_id[ii] == this_mpi_proc)
        {
            nodelist_owned[kk++] = ii;
        }
    }

    MPI_Barrier(MPI_COMM_WORLD);

    //vector<int>  nNode_owned_vector(n_mpi_procs), nNode_owned_sum(n_mpi_procs);
    vector<int>  nNode_owned_sum(n_mpi_procs);
    nNode_owned_vector.resize(n_mpi_procs);
    MPI_Allgather(&nNode_owned, 1, MPI_INT, &nNode_owned_vector[0], 1, MPI_INT, MPI_COMM_WORLD);

    nNode_owned_sum = nNode_owned_vector;
    for(int ii=1; ii < n_mpi_procs; ii++)
    {
        nNode_owned_sum[ii] += nNode_owned_sum[ii-1];
    }

    node_start = 0;
    node_end = 0;
    if(this_mpi_proc > 0) {
        node_start = nNode_owned_sum[this_mpi_proc-1];
    }
    node_end = nNode_owned_sum[this_mpi_proc]-1;

    LOG_INFO("pid: "      + to_string(this_mpi_proc) +
        "\tnNode_owned: " + to_string(nNode_owned) +
        "\tnode_start: "  + to_string(node_start) +
        "\tnode_end: "    + to_string(node_end + 1) +
        string(__FILE__) + ":" + to_string(__LINE__));

    MPI_Barrier(MPI_COMM_WORLD);

    //std::vector<int>  displs(n_mpi_procs);
    displs.resize(n_mpi_procs);
    displs[0] = 0;
    for(int ii=0; ii < n_mpi_procs-1; ii++) {
        displs[ii+1] = displs[ii] + nNode_owned_vector[ii];
    }

    MPI_Allgatherv(&nodelist_owned[0], nNode_owned, MPI_INT, &node_map_get_old[0], &nNode_owned_vector[0], &displs[0], MPI_INT, MPI_COMM_WORLD);


    for(int ii = 0; ii < nNode_global; ii++) {
        int n1 = node_map_get_old[ii];
        node_map_get_new[n1] = ii;
    }

    for(int ee = 0; ee < nElem_local; ee++) {
        int npElem = elemConn[ee].size()-3;

        for(int ii = 0; ii < npElem; ii++) {
            // LOG_INFO("elemConn: "      + to_string(elemConn[ee][3+ii]) +
            // string(__FILE__) + ":" + to_string(__LINE__));
            elemConn[ee][3+ii] = node_map_get_new[elemConn[ee][3+ii]];
            // LOG_INFO("elemConn: "      + to_string(elemConn[ee][3+ii]) +
            // string(__FILE__) + ":" + to_string(__LINE__));
        }
    }

    for(int ii=0; ii < NodalBCs.size(); ii++)
    {
        NodalBCs[ii].nodeNum = node_map_get_new[NodalBCs[ii].nodeNum];
    }

    for(int ii = 0; ii < NodalForces.size(); ii++)
    {
        NodalForces[ii].nodeNum = node_map_get_new[NodalForces[ii].nodeNum];
    }

    GeomData.node_map_get_old = node_map_get_old;
    GeomData.node_map_get_new = node_map_get_new;
    //GeomData.updateNodesAfterPartition();

    MPI_Barrier(MPI_COMM_WORLD);

    err_code = COMMON_SUCCESS;

    return err_code;
}


int  femSolids::readSurfaceElements(ifstream& infile, string& line)
{
    // read the number of elements
    getline(infile,line);    boost::trim(line);
/*
    nElemFaces = stoi(line);

    cout << "nElemFaces = " << nElemFaces << endl;

    vector<string>  stringlist;

    elemConnFaces.resize(nElemFaces);

    for(int ee=0; ee<nElemFaces; ee++)
    {
        getline(infile,line);
        boost::trim(line);
        //cout << line << endl;
        boost::algorithm::split(stringlist, line, boost::is_any_of("\t "), boost::token_compress_on);
        //stringlist.erase(std::remove(stringlist.begin(), stringlist.end(), " "), stringlist.end());


        int ind = stringlist.size();
        //cout << "ind = " << ind << endl;

        elemConnFaces[ee].resize(ind-1);

        for(int j=1; j<ind; j++)
            elemConnFaces[ee][j-1] = stoi(stringlist[j]) - 1;

        //printVector(elemConnFaces[ee]);
    }
*/
    return 0;
}



int  femSolids::readPrescribedBCs(ifstream& infile, string& line)
{
    // read the number of prescribed BCs
    getline(infile,line);    boost::trim(line);

    int nDBC = stoi(line), count;
    vector<string>  stringlist;

    cout << "nDBC = " << nDBC << endl;

    NodalBCs.resize(nDBC);
    for(int ee=0; ee<nDBC; ee++)
    {
        getline(infile,line);
        boost::trim(line);
        boost::algorithm::split(stringlist, line, boost::is_any_of("\t "), boost::token_compress_on);
        for(auto& str: stringlist)  boost::trim(str);

        count = stringlist.size();
        if(count < 3)
        {
            cerr << "Input error in 'prescribed boundary conditions' !" << endl;
        }


        NodalBCs[ee].nodeNum      = stoi(stringlist[0])-1 ;
        NodalBCs[ee].dof          = stoi(stringlist[1])-1 ;
        NodalBCs[ee].value        = stod(stringlist[2]) ;
        NodalBCs[ee].timeFuncNum  = -1;
        NodalBCs[ee].spaceFuncNum = -1;

        if(count >= 4)
          NodalBCs[ee].timeFuncNum  = stoi(stringlist[3])-1 ;

        if(count >= 5)
          NodalBCs[ee].spaceFuncNum = stoi(stringlist[4])-1 ;
    }

    return 0;
}





int  femSolids::readNodalForces(ifstream& infile, string& line)
{
    // read the number of nodal forces
    getline(infile,line);    boost::trim(line);

    int  nFBC = stoi(line), count;
    vector<string>  stringlist;

    cout << "nFBC = " << nFBC << endl;

    NodalForces.resize(nFBC);
    for(int ee=0; ee<nFBC; ee++)
    {
        getline(infile,line);
        boost::trim(line);
        boost::algorithm::split(stringlist, line, boost::is_any_of("\t "), boost::token_compress_on);
        for(auto& str: stringlist)  boost::trim(str);

        count = stringlist.size();
        if(count < 3)
        {
            cerr << "Input error in 'nodal forces' !" << endl;
        }


        NodalForces[ee].nodeNum      = stoi(stringlist[0])-1 ;
        NodalForces[ee].dof          = stoi(stringlist[1])-1 ;
        NodalForces[ee].value        = stod(stringlist[2]) ;
        NodalForces[ee].timeFuncNum  = -1;

        if(count >= 4)
          NodalForces[ee].timeFuncNum  = stoi(stringlist[3])-1 ;
    }

    return 0;
}






int femSolids::readMaterialProps(ifstream& infile, string& line)
{

    // read {
    getline(infile,line);    boost::trim(line);

    // read id
    getline(infile,line);    boost::trim(line);

    vector<string>  stringlist;
    boost::algorithm::split(stringlist, line, boost::is_any_of(":"), boost::token_compress_on);
    for(auto& str: stringlist)  boost::trim(str);

    int  id = stoi(stringlist[1]);

    // read name
    getline(infile,line);    boost::trim(line);
    boost::algorithm::split(stringlist, line, boost::is_any_of(":"), boost::token_compress_on);
    for(auto& str: stringlist)  boost::trim(str);


    std::string matkey(stringlist[1]);

    int  idd = getMaterialID(matkey);

    cout << " idd = " << idd << endl;

    if( (idd == 101) || (idd == 102) || (idd == 103) || (idd == 2006) || (idd == 2007) )
      intVarFlag = true;

    MatlDataList.push_back( NewMaterial(idd) );

    MatlDataList[numMaterials]->setID(id);

    MatlDataList[numMaterials]->SolnData = &(SolnData);

    MatlDataList[numMaterials]->readData(infile, line);

    numMaterials++;

    cout << "numMaterials = " << numMaterials << endl;

    return 0;
}





int femSolids::readElementProps(ifstream& infile, string& line)
{
    cout <<  " femSolids::readElementProps ... " << endl;

    ElementTypeDataList.push_back( new ElementTypeData );

    ElementTypeDataList[numElemTypes]->readData(infile, line);

    numElemTypes++;

    cout << "numElemTypes = " << numElemTypes << endl;

    cout <<  " femSolids::readElementProps ... " << endl;

    return 0;
}





int  femSolids::readBodyForce(ifstream& infile, string& line)
{
    vector<string>  stringlist;
    myPoint  bforce;

    // read {
    getline(infile,line);    boost::trim(line);

    while( infile && (line != "}") )
    {
        getline(infile,line);    boost::trim(line);

        if(line[0] == '}') break;

        if( isAValidLine(line) )
        {
            boost::algorithm::split(stringlist, line, boost::is_any_of(":"), boost::token_compress_on);
            for(auto& str: stringlist)  boost::trim(str);

            if(stringlist[0] == "value")
            {
                boost::algorithm::split(stringlist, stringlist[1], boost::is_any_of("\t "), boost::token_compress_on);
                for(auto& str: stringlist)  boost::trim(str);

                bforce[0] = stod(stringlist[0]);
                bforce[1] = stod(stringlist[1]);
                bforce[2] = stod(stringlist[2]);

                GeomData.setBodyForce(bforce);
            }
            else if( (stringlist[0] == "timefunction") || (stringlist[0] == "TimeFunction") )
            {
                GeomData.setBodyForceTimeFunction( stoi(stringlist[1])-1 );
            }
            else
            {
                throw runtime_error("Option not available in femSolids::readBodyForce");
            }
        }
    }


    return 0;
}



int  femSolids::readSolverDetails(ifstream& infile, string& line)
{
    vector<string>  stringlist;

    // read {
    getline(infile,line);    boost::trim(line);

    while( infile && (line != "}") )
    {
        getline(infile,line);    boost::trim(line);

        if(line[0] == '}') break;


        if( isAValidLine(line) )
        {
            boost::algorithm::split(stringlist, line, boost::is_any_of(":"), boost::token_compress_on);
            for(auto& str: stringlist)  boost::trim(str);

            if(stringlist[0] == "solvertype")
            {
                if(stringlist[1] == "arclength")
                {
                  NONLINEAR_SOLVER_TYPE = SOLVER_TYPE_ARCLENGTH;
                  ARC_LENGTH = true;
                }
                else
                  NONLINEAR_SOLVER_TYPE = SOLVER_TYPE_NEWTONRAPHSON;
            }
            else if(stringlist[1] == "formulation")
            {
            }
            else if(stringlist[0] == "timescheme")
            {
                SolnData.timeIntegrationScheme = stringlist[1];
            }
            else if(stringlist[0] == "spectralRadius")
            {
                SolnData.spectralRadius = stod(stringlist[1]);
            }
            else if(stringlist[0] == "finalTime")
            {
                timeFinal = stod(stringlist[1]);
            }
            else if(stringlist[0] == "timeStep")
            {
                boost::algorithm::split(stringlist, stringlist[1], boost::is_any_of("\t "), boost::token_compress_on);
                for(auto& str: stringlist)  boost::trim(str);

                if(stringlist.size() == 1)
                {
                    myTime.set( stod(stringlist[0]), stod(stringlist[0]), stod(stringlist[0]) );
                }
                else if(stringlist.size() == 2)
                {
                    myTime.set( stod(stringlist[0]), stod(stringlist[1]), stod(stringlist[0]) );
                }
                else
                {
                    myTime.set( stod(stringlist[0]), stod(stringlist[1]), stod(stringlist[2]) );
                }
            }
            else if(stringlist[0] == "maximumSteps")
            {
                stepsMax = stoi(stringlist[1]);
            }
            else if(stringlist[0] == "maximumIterations")
            {
                iterationsMax = stoi(stringlist[1]);
            }
            else if(stringlist[0] == "tolerance")
            {
                conv_tol = stod(stringlist[1]);
            }
            else if(stringlist[0] == "debug")
            {
                debug = ( stoi(stringlist[1]) == 1);
            }
            else if(stringlist[0] == "outputFrequency")
            {
                outputFreq = stoi(stringlist[1]);
            }
            else
            {
                throw runtime_error("Option not available in femSolids::readSolverDetails");
            }
        }
    }


    return 0;
}



int  femSolids::readTimeFunctions(ifstream& infile, string& line)
{
    // read {
    //getline(infile,line);    boost::trim(line);
    //prgReadTimeFunctions(infile);

    vector<string>  stringlist;
    vector<double>  doublelist;

    // read {
    getline(infile,line);    boost::trim(line);

    while( infile && (line != "}") )
    {
        getline(infile,line);    boost::trim(line);

        cout << line << endl;

        if(line[0] == '}') break;


        if( isAValidLine(line) )
        {
            boost::algorithm::split(stringlist, line, boost::is_any_of("\t "), boost::token_compress_on);
            for(auto& str: stringlist)  boost::trim(str);

            int id = stoi(stringlist[0]);

            if(id == 0)
            {
                cout << " Time function id '0' is not addmissible. It should be 1 or higher. " << endl;
                throw runtime_error(" Error in femSolids::readTimeFunctions...  ");
            }

            if(stringlist.size() < 11)
            {
                cout << " Number of temrs in Time function should be at least 11 " << endl;
                throw runtime_error(" Error in femSolids::readTimeFunctions...  ");
            }

            doublelist.resize(stringlist.size()-1);
            for(int i=0; i<stringlist.size()-1; ++i)
            {
                doublelist[i] = stod(stringlist[i+1]);
            }

            id -= 1;
            if(id == timeFunction.size())
            {
                timeFunction.push_back(make_unique<TimeFunction>());
                timeFunction[id]->setID(id);
            }
            timeFunction[id]->addTimeBlock(doublelist);
        }
    }

    for(auto& tmf : timeFunction)
    {
        tmf->printSelf();
        tmf->update();
    }

    return 0;
}





int femSolids::readInitialConditions(ifstream& infile, string& line)
{
    vector<string>  stringlist;
    string label;

    return 0;
}





int femSolids::readNodalDataForOutput(ifstream& infile, string& line)
{
/*
    // read the number of rows
    getline(infile,line);    boost::trim(line);

    int  nOut = stoi(line);

    vector<string>  stringlist;
    vector<double>  vecTemp(3);

    for(int ee=0; ee<nOut; ee++)
    {
        getline(infile,line);        boost::trim(line);
        boost::algorithm::split(stringlist, line, boost::is_any_of("\t "), boost::token_compress_on);
        //stringlist.erase(std::remove(stringlist.begin(), stringlist.end(), " "), stringlist.end());

        int ind = stringlist.size();
        //cout << "ind = " << ind << endl;

        vecTemp[0] = stod(stringlist[0]);
        vecTemp[1] = stod(stringlist[1]);
        vecTemp[2] = stod(stringlist[2]);

        OutputData.push_back(vecTemp);

        //printVector(DirichletBCs[ee]);
    }
*/
    return 0;
}


int femSolids::readOutputDetails(ifstream& infile, string& line)
{
    vector<string>  stringlist;

    // read {
    getline(infile,line);    boost::trim(line);

    while( infile && (line != "}") )
    {
        getline(infile,line);    boost::trim(line);

        cout << line << endl;

        if(line[0] == '}') break;


        if( isAValidLine(line) )
        {
            boost::algorithm::split(stringlist, line, boost::is_any_of(":"), boost::token_compress_on);
            for(auto& str: stringlist)  boost::trim(str);

            if(stringlist[0] == "Frequency")
            {
                outputfrequency = stoi(stringlist[1]);
            }
            else if(stringlist[0] == "Nodal")
            {
                outputlist_nodal.push_back(stringlist[1]);
            }
            else if(stringlist[0] == "Elemental")
            {
                outputlist_elemental.push_back(stringlist[1]);
            }
            else
            {
                throw runtime_error("Option not available in femSolids::readOutputDetails");
            }
        }
    }

    return 0;
}




int femSolids::printInfo()
{
    return 0;
}







int femSolids::setSpecifiedDOFs(vector<vector<bool> >& NodeType)
{
    int  ii, nn, dof;

    dofs_specified.clear();
    for(ii=0; ii<NodalBCs.size(); ii++)
    {
        nn = NodalBCs[ii].nodeNum;
        dof = NodalBCs[ii].dof;

        NodeType[nn][dof] = true;
        dofs_specified.push_back(nn*ndof+dof);
    }
    findUnique(dofs_specified);

    return 0;
}







int femSolids::setBoundaryConditions()
{
    int  ii, nn, dof, timeFuncNum;
    double loadFactor;

    //SolnData.dispApplied.setZero();

    //wj

    SolnData.dispAppliedLocal.setZero();
    for(ii=0; ii<NodalBCs.size(); ii++)
    {
        nn          = NodalBCs[ii].nodeNum;
        dof         = NodalBCs[ii].dof;
        timeFuncNum = NodalBCs[ii].timeFuncNum;

        loadFactor = 1.0;
        if(timeFuncNum >= 0)
          loadFactor = timeFunction[timeFuncNum]->getFactor();

        //SolnData.dispApplied[nn*ndof+dof] = NodalBCs[ii].value * loadFactor;
        // LOG_INFO("SolnData.dispApplied"      + to_string(SolnData.dispApplied[nn*ndof+dof]) + ":" + to_string(nn) + 
        // string(__FILE__) + ":" + to_string(__LINE__));

        if(node_map_g2l.find(nn) != node_map_g2l.end())
        {
            SolnData.dispAppliedLocal[node_map_g2l[nn] * ndof + dof] = NodalBCs[ii].value * loadFactor;
            // LOG_INFO("SolnData.dispAppliedLocal"      + to_string(SolnData.dispAppliedLocal[node_map_g2l[nn] * ndof + dof]) + ":" + to_string(nn) + 
            // string(__FILE__) + ":" + to_string(__LINE__));
        }

    }

      
   // SolnData.dispApplied -= SolnData.disp;

    SolnData.dispAppliedLocal -= SolnData.dispLocal;


/*
    // spatially varying boundary conditions
    nn = nodeNums[i];

    xc = GeomData.NodePosOrig[nn][0];
    yc = GeomData.NodePosOrig[nn][1];
    zc = GeomData.NodePosOrig[nn][2];

    value = mathfun.getValue(xc, yc, zc) * loadFactor;
*/

    return 0;
}




int femSolids::addBoundaryConditions()
{
    int ii, ind;
    for(ii=0; ii<NodalBCs.size(); ++ii)
    {
        ind = NodalBCs[ii].nodeNum * ndof + NodalBCs[ii].dof;

        SolnData.disp[ind]  += SolnData.dispApplied[ind];
    }

    return 0;
}







int femSolids::setInitialConditions()
{
    PetscPrintf(MPI_COMM_WORLD, "femSolids::setInitialConditions() ... \n");
/*
    double  x=0.0, y=0.0, z=0.0, eps = 0.02;
    string  expr;

    SolnData.dispPrev.setZero();
    
    for(int dof=0; dof<ndof; ++dof)
    {
      expr = initialConitions[dof];

      cout << expr << endl;

      if( !expr.empty() )
      {
        myMathFunction  mathfun;
        mathfun.initialise(initialConitions[dof]);

        exprtk::expression<double>   expression;

        exprtk::symbol_table<double> symbol_table;
        exprtk::parser<double>       parser;

        symbol_table.add_variable("x",x);
        symbol_table.add_variable("y",y);
        symbol_table.add_variable("z",z);
        symbol_table.add_variable("eps",eps);
        symbol_table.add_constants();
        expression.register_symbol_table(symbol_table);

        parser.compile(expr, expression);

        for(int ii=0; ii<nNode_global; ++ii)
        {
          x = nodeCoordsOrig[ii][0];
          y = nodeCoordsOrig[ii][1];
          z = nodeCoordsOrig[ii][2];

          SolnData.dispPrev(node_map_get_new[ii]*ndof + dof) = expression.value();
        }
      }
    }
    PetscPrintf(MPI_COMM_WORLD, "femSolids::setInitialConditions() ... \n");

    // add boundary Conditions
    //setBoundaryConditions();
    //solnPrev += solnApplied;

    SolnData.disp = SolnData.dispPrev;
*/

    SolnData.velo.setZero();
    SolnData.veloPrev = SolnData.velo;

    PetscPrintf(MPI_COMM_WORLD, "femSolids::setInitialConditions() ... \n");

/*
    double  xx=0.0, yy=0.0, zz=0.0, fact;

    for(int ii=0; ii<nNode; ++ii)
    {
        xx = nodeCoordsOrig[ii][0];
        yy = nodeCoordsOrig[ii][1];
        //zz = nodeCoordsOrig[ii][2];

        //veloPrev(ii*2) = 2.0*yy*(3.0-yy)/3.0;
        SolnData.veloPrev(ii*ndim) = 1.0*yy;

        //veloPrev(ii*ndim) = 16.0*0.45*yy*zz*(0.41-yy)*(0.41-zz)/0.41/0.41/0.41/0.41;
    }
    SolnData.velo = SolnData.veloPrev;
*/

    return 0;
}




int femSolids::addExternalForces()
{
    solverPetsc->Fext.setZero();

    if(this_mpi_proc > 0) return 0;

    if(debug) {PetscPrintf(MPI_COMM_WORLD, "\n\n femSolids::addExternalForces() ... STARTED \n");}

    if(NodalForces.size() > 0)
    {
        int  nn, dof, ii, timeFuncNum;

        for(ii=0;ii<NodalForces.size();ii++)
        {
            nn          = NodalForces[ii].nodeNum;
            dof         = NodalForces[ii].dof;
            timeFuncNum = NodalForces[ii].timeFuncNum;

            loadFactor = 1.0;
            if(timeFuncNum >= 0)
              loadFactor = timeFunction[timeFuncNum]->getFactor();

            solverPetsc->Fext[nn*ndof+dof] += loadFactor*NodalForces[ii].value;
        }
    }

    if(debug) {PetscPrintf(MPI_COMM_WORLD, "\n\n\n femSolids::addExternalForces() ... ENDED \n");}

    return 0;
}






int femSolids::timeUpdate()
{
    if(debug) {PetscPrintf(MPI_COMM_WORLD, " femSolids::timeUpdate() ... STARTED \n");}

    myTime.update();

    // set parameters for the time integration scheme.
    // need to be done every time step to account for adaptive time stepping
    SolnData.setTimeParam();

    int  idd = 0;//SolnData.ElemProp[0]->id;
    if( (idd == 2010) || (idd == 2060) )
    {
        for(int ee=0; ee<nElem_local; ee++)
        {
            elems[ee]->presDOF     = elems[ee]->presDOFprev;
            elems[ee]->presDOFprev = elems[ee]->presDOFprev2;
        }
    }

    // copy internal variables intVar1 to intVar2
    if(intVarFlag)
      copyElemInternalVariables();

    // update time functions
    for(auto& tmf : timeFunction)
      tmf->update();

    // set Dirichlet boundary conditions
    setBoundaryConditions();

    if(debug) {PetscPrintf(MPI_COMM_WORLD, "\n femSolids::timeUpdate() ... FINISHED \n\n"); }

    return 0;
}




int femSolids::updateIterStep()
{
    SolnData.updateIterStep();

    int  ii, jj, ind;
    // for(ii=0;ii<nNode_global;ii++)
    // {
    //   for(jj=0;jj<ndim;jj++)
    //   {
    //     ind = ii*ndof+jj;
    //     GeomData.NodePosNew[ii][jj] = GeomData.NodePosOrig[ii][jj] + SolnData.disp[ind];
    //     // LOG_INFO("GeomData.NodePosNew1 = " + to_string(GeomData.NodePosOrig[ii][jj]) +  ":"+ to_string(SolnData.disp[ind]) +  ":"+
    //     // to_string(ii) + string(__FILE__) + ":" + to_string(__LINE__));
    //     GeomData.NodePosCur[ii][jj] = GeomData.NodePosOrig[ii][jj] + SolnData.dispCur[ind];
    //     //LOG_INFO("GeomData.NodePosNew = " + to_string(GeomData.NodePosNew[ii][jj]) +  ":"+ to_string(GeomData.NodePosCur[ii][jj]) +  ":"+
    //     //to_string(ii) + string(__FILE__) + ":" + to_string(__LINE__));
    //   }
    // }


    for(ii=0;ii<num_nodeInLocalElem;ii++)
    {
      for(jj=0;jj<ndim;jj++)
      {
        //ind = node_map_l2g[ii] * ndof + jj;

        ind = ii * ndof + jj;
        GeomData.NodePosNewLocal[ii][jj] = GeomData.NodePosOrigLocal[ii][jj] + SolnData.dispLocal[ind];

        // LOG_INFO("GeomData.NodePosNewLocal =" + to_string(GeomData.NodePosOrigLocal[ii][jj]) +  ":"+ to_string(SolnData.disp[ind]) +  ":"+ 
        // to_string(node_map_l2g[ii]) + string(__FILE__) + ":" + to_string(__LINE__));
        // LOG_INFO("GeomData.NodePosNewLocal =" + to_string(GeomData.NodePosNewLocal[ii][jj]) +  ":" + to_string(GeomData.NodePosOrigLocal[ii][jj]) +  ":"+ to_string(SolnData.disp[ind]) +  ":"+ to_string(ind) +  ":"+ to_string(ii * ndof + jj) +  ":" +  to_string(SolnData.dispLocal[ii * ndof + jj])+
        // to_string(node_map_l2g[ii]) + string(__FILE__) + ":" + to_string(__LINE__));

        GeomData.NodePosCurLocal[ii][jj] = GeomData.NodePosOrigLocal[ii][jj] + SolnData.dispCurLocal[ind];
        //LOG_INFO("GeomData.NodePosNewLocal =" + to_string(GeomData.NodePosNewLocal[ii][jj]) +  ":"+ to_string(GeomData.NodePosCurLocal[ii][jj]) +  ":"+ 
        //to_string(node_map_l2g[ii]) + string(__FILE__) + ":" + to_string(__LINE__));

        // LOG_INFO("GeomData.NodePosCurLocal =" + to_string(GeomData.NodePosCurLocal[ii][jj]) +  ":" + to_string(GeomData.NodePosOrigLocal[ii][jj]) +  ":"+ to_string(SolnData.dispCur[ind]) +  ":"+  to_string(SolnData.dispCurLocal[ii * ndof + jj])+
        // to_string(node_map_l2g[ii]) + string(__FILE__) + ":" + to_string(__LINE__));

      }
    }

    return 0;
}




int  femSolids::reset()
{
    SolnData.reset();

    return 0;
}




bool femSolids::converged()
{
  if (rhsNorm < conv_tol)
    return true;

  return false;
}



bool femSolids::diverging(double factor)
{
  if (isnan(rhsNorm)) return true;

  if (rhsNormPrev > -0.1 && (rhsNorm/rhsNormPrev) > factor) return true;

  if (localStiffnessError != 0) return true;

  return false;
}




int  femSolids::saveSolution()
{
    SolnData.saveSolution();

    int  idd = 0.0;//SolnData.ElemProp[0]->id;
    if( (idd == 2010) || (idd == 2060) )
    {
        for(int ee=0; ee<nElem_local; ee++)
        {
            elems[ee]->presDOFprev2 = elems[ee]->presDOFprev;
            elems[ee]->presDOFprev  = elems[ee]->presDOF;
        }
    }

    if(intVarFlag)
    {
        for(int ee=0; ee<nElem_local; ee++)
        {
          elems[ee]->ivar.saveSolution();
        }
    }

    return 0;
}




int  femSolids::copyElemInternalVariables()
{
    return 0;
}





int femSolids::writeNodalData()
{
  return 0;
}





int  femSolids::readInputBDF(ifstream& infile)
{
    cout << " femSolids::readInputBDF " << endl;

    string line, line2;
    vector<string>  stringlist, stringlist2;

    vector<vector<int> >  elemConnTemp;
    vector<myPoint>       nodecoords;
    vector<int>  vecTempInt;
    vector<double>   vecTempDbl;
    myPoint   pt1, pt2;
    
    int  ee, ii, jj, kk;
    
    nNode_global = 0;
    nElem_global = 0;
    ndim = 3;
    ndof = 3;


    while(getline(infile,line))
    {
      //std::cout << line << std::endl;

      boost::algorithm::split(stringlist, line, boost::is_any_of(" "), boost::token_compress_on);
      for(auto& str: stringlist)  boost::trim(str);


      if( (stringlist[0] == "GRID") || (stringlist[0] == "GRID*") )
      {
          if(stringlist[0] == "GRID")
          {
            pt1[0] = stod(stringlist[2]);
            pt1[1] = stod(stringlist[3]);
            pt1[2] = stod(stringlist[4]);
          }
          else
          {
            // 0-7, 8-23, 24-39, 40-55, 56-72

            pt1[0] = stod( line.substr(40,16) );
            pt1[1] = stod( line.substr(56,16) );

            getline(infile,line2);
            pt1[2] = stod( line2.substr(8,16) );
          }

          nodecoords.push_back(pt1);
      }
      //
      //
      else if(stringlist[0] == "CQUAD4")
      {
        ndim = 2;
        ndof = 2;
        //cout << " size = " << stringlist.size() << endl;

        vecTempInt.clear();

        for(ii=0; ii<stringlist.size()-1; ii++)
        {
            vecTempInt.push_back( stoi(stringlist[1+ii]) );
        }

        elemConn.push_back(vecTempInt);
      }
      //
      //
      else if(stringlist[0] == "MAT1") // Isotropic linear elastic material
      {
          MatlDataList.push_back( NewMaterial(1) );

          MatlDataList[numMaterials]->setID( stoi(stringlist[1]) );

          MatlDataList[numMaterials]->SolnData = &(SolnData);
          
          vecTempDbl.clear();
          for(ii=0; ii<stringlist.size()-2; ii++)
              vecTempDbl.push_back( stod(stringlist[2+ii]) );

          MatlDataList[numMaterials]->setData(vecTempDbl);

          numMaterials++;

          cout << "numMaterials = " << numMaterials << endl;
      }
      //
      //
      else if(stringlist[0] == "SPC1")
      {
      }
    }
    
    nNode_global = nodecoords.size();
    nElem_global = elemConn.size();

    GeomData.setDimension(ndim);
    GeomData.setNodalPositions(nodecoords);


    return 0;
}
























